$Re$はレイノルズ数.
(2)式の発散をとって, \begin{equation} \frac{\partial D}{\partial t} + \mathrm{div} ((\boldsymbol{v}\cdot\nabla)\boldsymbol{v}) = - \Delta p + \frac{1}{Re}\Delta D, \quad D \equiv \nabla \cdot \boldsymbol{v} \notag \end{equation} \begin{equation} \therefore \: \Delta p = -\frac{\partial D}{\partial t} - \mathrm{div} ((\boldsymbol{v}\cdot\nabla)\boldsymbol{v}) + \frac{1}{Re}\Delta D \end{equation} $D=0$という条件を使えば, \begin{equation} \therefore \: \Delta p = - \mathrm{div} ((\boldsymbol{v}\cdot\nabla)\boldsymbol{v}) \tag{3} \end{equation} というポアソン方程式を得る.しかし数値的には必ずしも$D=0$になるとは限らないから時間微分項は,数値計算するときを想定して,単純に
\begin{equation} \frac{\partial D}{\partial t} \simeq \frac{D_{n+1} - D_n}{\Delta t} \end{equation}としておいて,次の時刻$t_{n+1}$で$D$がゼロになるように
\begin{equation} \frac{\partial D}{\partial t} \simeq -\frac{D_n}{\Delta t} \end{equation}というふうにする.あと実際$D$は小さい数になると思われるので2階微分はほぼゼロとして消す.なお,
\begin{align} \Delta\boldsymbol{v} &=\mathrm{grad}\:\mathrm{div}\boldsymbol{v}-\mathrm{rot}\:\mathrm{rot}\boldsymbol{v} \\ \therefore \: \mathrm{div} \Delta\boldsymbol{v} &= \mathrm{div} \mathrm{grad}(\mathrm{div}\boldsymbol{v})-\mathrm{div} \mathrm{rot} \:(\mathrm{rot}\boldsymbol{v}) \\ &= \Delta (\mathrm{div}\boldsymbol{v}) = \Delta D \end{align}を用いた.
(円柱に近いほど粗くなる格子を設定している.)
まず普通の極座標変換$ x = r \cos\theta, \: y = r \sin\theta$をやってから,$ r = \mathrm{e}^{\xi}$, $\displaystyle \frac{\partial}{\partial r} = \frac{\partial \xi}{\partial r} \frac{\partial}{\partial \xi} = \mathrm{e}^{-\xi} \frac{\partial}{\partial \xi}$の変換を行えばいい.
半径方向の速度を$v_r$,動径方向の速度を$v_\theta$とする.すなわち,$ \boldsymbol{v} = v_r \vec{e}_r + v_\theta \vec{e}_\theta$としておく.
と書けるが,$\vec{e}_r,\vec{e}_\theta$は$\theta$によって変化するので,この単位ベクトルの動径方向微分を考えなければならない.
\begin{align} \vec{e}_r &= \cos\theta\vec{e}_x + \sin\theta\vec{e}_y \\ \vec{e}_\theta &= -\sin\theta\vec{e}_x + \cos\theta\vec{e}_y \end{align}より,
\begin{align} \frac{\partial}{\partial \theta} \vec{e}_r &= -\sin\theta\vec{e}_x + \cos\theta\vec{e}_y = \vec{e}_\theta \\ \frac{\partial}{\partial \theta} \vec{e}_\theta &= -\cos\theta\vec{e}_x + \sin\theta\vec{e}_y = -\vec{e}_r \end{align}だから,これより
\begin{align} (\vec{v}\cdot\nabla)\vec{v} &= \left( v_r\frac{\partial}{\partial r} + \frac{v_\theta}{r} \frac{\partial}{\partial \theta} \right) (v_r \vec{e}_r + v_\theta \vec{e}_\theta) \\ &=v_r\frac{\partial v_r}{\partial r}\vec{e}_r + \frac{v_\theta}{r} \frac{\partial v_r }{\partial \theta}\vec{e}_r + \frac{v_r v_\theta}{r} \frac{\partial \vec{e}_r }{\partial \theta} \\ &\quad + v_r\frac{\partial v_\theta}{\partial r}\vec{e}_\theta + \frac{v_\theta}{r} \frac{\partial v_\theta }{\partial \theta}\vec{e}_\theta + \frac{v_\theta^2}{r} \frac{\partial \vec{e}_\theta }{\partial \theta} \\ &= \vec{e}_r \left( v_r\frac{\partial v_r}{\partial r} + \frac{v_\theta}{r} \frac{\partial v_r }{\partial \theta} - \frac{v_\theta^2}{r} \right) + \vec{e}_\theta \left( v_r\frac{\partial v_\theta}{\partial r} + \frac{v_\theta}{r} \frac{\partial v_\theta }{\partial \theta} + \frac{v_r v_\theta}{r} \right) \\ &= \mathrm{e}^{-\xi}\left( v_r\frac{\partial v_r}{\partial \xi} + v_\theta \frac{\partial v_r }{\partial \theta} - v_\theta^2 \right)\vec{e}_r + \mathrm{e}^{-\xi} \left( v_r\frac{\partial v_\theta}{\partial \xi} + v_\theta\frac{\partial v_\theta }{\partial \theta} + v_r v_\theta \right)\vec{e}_\theta \end{align}と書ける.(移流項の計算は解析力学のEuler-Lanrange方程式を使って求めることもできる.) ということで,これでやっと(2)式が変換できる.
\begin{align} \frac{\partial v_r}{\partial t} + \mathrm{e}^{-\xi}\left( v_r\frac{\partial v_r}{\partial \xi} + v_\theta \frac{\partial v_r }{\partial \theta} - v_\theta^2 \right) &= -\mathrm{e}^{-\xi} \frac{\partial p}{\partial \xi} + \frac{\mathrm{e}^{-2\xi}}{Re} \left( \frac{\partial^2 v_r}{\partial \xi^2} + \frac{\partial^2 v_r}{\partial \theta^2} - v_r- 2\frac{\partial v_\theta}{\partial \theta} \right) \\ \frac{\partial v_\theta}{\partial t} + \mathrm{e}^{-\xi} \left( v_r\frac{\partial v_\theta}{\partial \xi} + v_\theta\frac{\partial v_\theta }{\partial \theta} + v_r v_\theta \right) &= -\mathrm{e}^{-\xi} \frac{\partial p}{\partial \theta} + \frac{\mathrm{e}^{-2\xi}}{Re} \left( \frac{\partial^2 v_\theta}{\partial \xi^2} + \frac{\partial^2 v_\theta}{\partial \theta^2} - v_\theta+ 2\frac{\partial v_r}{\partial \theta} \right) \end{align}連続の式:
\begin{equation} \nabla\cdot \vec{v} = 0 \Leftrightarrow \frac{1}{r} \frac{\partial}{\partial r} (r v_r) + \frac{1}{r} \frac{\partial v_\theta}{\partial \theta} = \frac{\partial v_r}{\partial r} + \frac{v_r}{r} + \frac{1}{r} \frac{\partial v_\theta}{\partial \theta} = 0 \end{equation}を用いて $(\star), (\star\star)$の部分を変換する.
\begin{align} (\star) &= \frac{v_r}{r} \frac{\partial}{\partial r} \left( v_r + \frac{\partial v_\theta}{\partial \theta} \right) + v_r \frac{\partial^2 v_r}{\partial r^2} + \frac{v_r}{r^2} \frac{\partial v_\theta}{\partial \theta} \\ &= \frac{v_r}{r} \frac{\partial}{\partial r} \left( -r \frac{\partial v_r}{\partial r} \right) + v_r \frac{\partial^2 v_r}{\partial r^2} + \frac{v_r}{r^2} \frac{\partial v_\theta}{\partial \theta} \\ &= -v_r \frac{\partial^2 v_r}{\partial r^2} - \frac{v_r}{r} \frac{\partial v_r}{\partial r} + v_r \frac{\partial^2 v_r}{\partial r^2} + \frac{v_r}{r^2} \frac{\partial v_\theta}{\partial \theta} \\ &= -\frac{v_r}{r} \frac{\partial v_r}{\partial r} + \frac{v_r}{r^2} \frac{\partial v_\theta}{\partial \theta} \\ (\star\star) &= \frac{v_\theta}{r^2} \frac{\partial}{\partial \theta} \left( v_r + \frac{\partial v_\theta}{\partial \theta} \right) + \frac{v_\theta}{r} \frac{\partial^2 v_r}{\partial r \partial \theta} - \frac{2 v_\theta}{r} \frac{\partial v_\theta}{\partial r} \\ &= \frac{v_\theta}{r^2} \frac{\partial}{\partial \theta} \left( -r \frac{\partial v_r}{\partial r} \right) + \frac{v_\theta}{r} \frac{\partial^2 v_r}{\partial r \partial \theta} - \frac{2 v_\theta}{r} \frac{\partial v_\theta}{\partial r} \\ &= -\frac{v_\theta}{r} \frac{\partial^2 v_r}{\partial r \partial\theta} + \frac{v_\theta}{r} \frac{\partial^2 v_r}{\partial r \partial \theta} - \frac{2 v_\theta}{r} \frac{\partial v_\theta}{\partial r} \\ &= -\frac{2 v_\theta}{r} \frac{\partial v_\theta}{\partial r} \end{align}より,
\begin{align} (\star) + (\star\star) &= -\frac{v_r}{r} \frac{\partial v_r}{\partial r} + \frac{v_r}{r^2} \frac{\partial v_\theta}{\partial \theta} -\frac{2 v_\theta}{r} \frac{\partial v_\theta}{\partial r} \\ &= \frac{v_r}{r} \left( \frac{v_r}{r} + \frac{1}{r}\frac{\partial v_\theta}{\partial \theta} \right) + \frac{v_r}{r^2} \frac{\partial v_\theta}{\partial \theta} -\frac{2 v_\theta}{r} \frac{\partial v_\theta}{\partial r} \\ &= \frac{v_r^2}{r^2} + \frac{2 v_r}{r^2} \frac{\partial v_\theta}{\partial \theta} -\frac{2 v_\theta}{r} \frac{\partial v_\theta}{\partial r} \end{align}が得られる.したがって,
\begin{align} \mathrm{div} ((\boldsymbol{v}\cdot\nabla)\boldsymbol{v}) &= \left( \frac{\partial v_r}{\partial r} \right)^2 + \frac{2}{r} \frac{\partial v_r}{\partial \theta} \frac{\partial v_\theta}{\partial r} + \frac{1}{r^2} \left(\frac{\partial v_\theta}{\partial \theta} \right)^2 \\ &\quad + \frac{v_r^2}{r^2} + \frac{2 v_r}{r^2} \frac{\partial v_\theta}{\partial \theta} -\frac{2 v_\theta}{r} \frac{\partial v_\theta}{\partial r} \\ & = \frac{1}{\mathrm{e}^{2\xi}} \left\{ \left( \frac{\partial v_r}{\partial \xi} \right)^2 + 2 \frac{\partial v_r}{\partial \theta} \frac{\partial v_\theta}{\partial \xi} + \left(\frac{\partial v_\theta}{\partial \theta} \right)^2 + v_r^2 + 2\left( v_r \frac{\partial v_\theta}{\partial \theta} - v_\theta \frac{\partial v_\theta}{\partial \xi} \right) \right\} \end{align}まで変形できた.以上より(3)式は
\begin{align} \frac{\partial^2 p}{\partial \xi^2} + \frac{\partial^2 p}{\partial \theta^2} = -\left\{ \left( \frac{\partial v_r}{\partial \xi} \right)^2 + 2 \frac{\partial v_r}{\partial \theta} \frac{\partial v_\theta}{\partial \xi} + \left(\frac{\partial v_\theta}{\partial \theta} \right)^2 + v_r^2 + 2\left( v_r \frac{\partial v_\theta}{\partial \theta} - v_\theta \frac{\partial v_\theta}{\partial \xi} \right) \right\} \end{align}というポアソン方程式に変換される。あとは数値計算のために$D$の時間微分の項を付け加えるだけ.
\begin{align} D = \nabla\cdot\vec{v} = \frac{\partial v_r}{\partial r} + \frac{v_r}{r} + \frac{1}{r} \frac{\partial v_\theta}{\partial \theta} = \mathrm{e}^{-\xi} \left( \frac{\partial v_r}{\partial \xi} + v_r + \frac{\partial v_\theta}{\partial \theta} \right) \end{align}なので,これに$\mathrm{e}^{2\xi}$をかけて$\Delta t$で割ったものを加えることで,数値計算に使う式は,
\begin{align} \frac{\partial^2 p}{\partial \xi^2} + \frac{\partial^2 p}{\partial \theta^2} = &-\left\{ \left( \frac{\partial v_r}{\partial \xi} \right)^2 + 2 \frac{\partial v_r}{\partial \theta} \frac{\partial v_\theta}{\partial \xi} + \left(\frac{\partial v_\theta}{\partial \theta} \right)^2 + v_r^2 + 2\left( v_r \frac{\partial v_\theta}{\partial \theta} - v_\theta \frac{\partial v_\theta}{\partial \xi} \right) \right\}\\ & + \frac{\mathrm{e}^{\xi} }{\Delta t} \left( \frac{\partial v_r}{\partial \xi} + v_r + \frac{\partial v_\theta}{\partial \theta} \right) \end{align}となる.
using Plots
Plots.pyplot()
Plots.PyPlotBackend()
i方向が$\xi$,j方向が$\theta$
function BCforP(p)
for j=1:JM
p[IM, j]=0.0
p[1,j]=(4*p[2,j]-p[3,j])/3.0
end
end
function BCforV(u,v, dtheta)
for j=1:JM
u[1,j]=0.0
v[1,j]=0.0
u[IM,j]=cos((j-1)*dtheta)
v[IM,j]=-sin((j-1)*dtheta)
end
end
BCforV (generic function with 1 method)
function SolvePoisson(p, u,v, dtheta, dxi, dt, rhs)
# rhs=zeros(IM,JM)
for i=2:IM-1
for j=2:JM-1
uxi=(u[i+1,j]-u[i-1,j])/(2*dxi)
vxi=(v[i+1,j]-v[i-1,j])/(2*dxi)
utheta=(u[i,j+1]-u[i,j-1])/(2*dtheta)
vtheta=(v[i,j+1]-v[i,j-1])/(2*dtheta)
rhs[i,j] = uxi^2 + 2*vxi*utheta + vtheta^2 + u[i,j]^2 + 2*(u[i,j]*vtheta-v[i,j]*vxi) -
exp((i-1)*dxi)/dt*(uxi+u[i,j]+vtheta)
end
end
for i=2:IM-1
uxi=(u[i+1,1]-u[i-1,1])/(2*dxi)
vxi=(v[i+1,1]-v[i-1,1])/(2*dxi)
utheta=(u[i,2]-u[i,JM-1])/(2*dtheta)
vtheta=(v[i,2]-v[i,JM-1])/(2*dtheta)
rhs[i,1] = uxi^2 + 2*vxi*utheta + vtheta^2 + u[i,1]^2 + 2*(u[i,1]*vtheta-v[i,1]*vxi) -
exp((i-1)*dxi)/dt*(uxi+u[i,1]+vtheta)
end
@.(rhs[2:IM-1,JM] = copy(rhs[2:IM-1,1]))
for iteration=1:MAXITP
res=0.0
for i=2:IM-1
dp= (dtheta^2*(p[i+1,1]+p[i-1,1]) + dxi^2*(p[i,2]+p[i,JM-1]) + dtheta^2*dxi^2*rhs[i,1])/2.0/(dxi^2+dtheta^2)
dp=dp-p[i,1]
res+=dp^2
p[i,1]+=OMEGA*dp
end
p[2:IM-1,JM]=copy(p[2:IM-1,1])
for i=2:IM-1
for j=2:JM-1
dp= (dtheta^2*(p[i+1,j]+p[i-1,j]) + dxi^2*(p[i,j+1]+p[i,j-1]) + dtheta^2*dxi^2*rhs[i,j])/2.0/(dxi^2+dtheta^2)
dp=dp-p[i,j]
res+=dp^2
p[i,j]+=OMEGA*dp
end
end
BCforP(p)
res=sqrt(res/(IM*JM))
if res < ERRORP
resp=res
itrp=iteration
break
end
end
end
SolvePoisson (generic function with 1 method)
function VelocityEQ(p,u,v,dtheta,dxi,dt,urhs,vrhs)
# urhs=zeros(IM,JM)
# vrhs=zeros(IM,JM)
for i=2:IM-1
urhs[i,1]= -exp(-(i-1)*dxi)*(p[i+1,1]-p[i-1,1])/(2*dxi) + exp(-2*(i-1)*dxi)/RE * ( (u[i+1,1]-2*u[i,1]+u[i-1,1])/(dxi^2) +
(u[i,2]-2*u[i,1]+u[i,JM-1])/(dtheta^2) - u[i,1] - (v[i,2]-v[i,JM-1])/dtheta )
vrhs[i,1]= -exp(-(i-1)*dxi)*(p[i,2]-p[i,JM-1])/(2*dtheta) + exp(-2*(i-1)*dxi)/RE *
( (v[i+1,1]-2*v[i,1]+v[i-1,1])/(dxi^2) + (v[i,2]-2*v[i,1]+v[i,JM-1])/(dtheta^2) - v[i,1] + (u[i,2]-u[i,JM-1])/dtheta )
urhs[i,JM]=urhs[i,1]
vrhs[i,JM]=vrhs[i,1]
end
for i=2:IM-1
for j=2:JM-1
urhs[i,j]= -exp(-(i-1)*dxi)*(p[i+1,j]-p[i-1,j])/(2*dxi) + exp(-2*(i-1)*dxi)/RE *
( (u[i+1,j]-2*u[i,j]+u[i-1,j])/(dxi^2) +
(u[i,j+1]-2*u[i,j]+u[i,j-1])/(dtheta^2) - u[i,j] - (v[i,j+1]-v[i,j-1])/dtheta )
vrhs[i,j]= -exp(-(i-1)*dxi)*(p[i,j+1]-p[i,j-1])/(2*dtheta) + exp(-2*(i-1)*dxi)/RE *
( (v[i+1,j]-2*v[i,j]+v[i-1,j])/(dxi^2) +
(v[i,j+1]-2*v[i,j]+v[i,j-1])/(dtheta^2) - v[i,j] + (u[i,j+1]-u[i,j-1])/dtheta )
end
end
for i=2:IM-1
for j=1:JM
urhs[i,j] -= exp(-(i-1)*dxi) * (-v[i,j]^2)
vrhs[i,j] -= exp(-(i-1)*dxi) * (v[i,j]*u[i,j])
end
end
#移流項 xi方向
for i=2:IM-1
for j=1:JM
if i==2
urhs[i,j] -= exp(-(i-1)*dxi) * ( u[i,j] * (-u[i+2,j]+8.0*(u[i+1,j]-u[i-1,j])+(2.0*u[i-1,j]-u[i,j]))/(12.0*dxi) +
abs(u[i, j]) * (u[i+2,j]-4.0*u[i+1,j]+6.0*u[i,j]-4.0*u[i-1,j]+(2.0*u[i-1,j]-u[i,j])) / (4.0*dxi))
vrhs[i,j] -= exp(-(i-1)*dxi) * ( u[i,j] * (-v[i+2,j]+8.0*(v[i+1,j]-v[i-1,j])+(2.0*v[i-1,j]-v[i,j]))/(12.0*dxi) +
abs(u[i, j]) * (v[i+2,j]-4.0*v[i+1,j]+6.0*v[i,j]-4.0*v[i-1,j]+(2.0*v[i-1,j]-v[i,j])) / (4.0*dxi))
elseif i==IM-1
urhs[i,j] -= exp(-(i-1)*dxi) * ( u[i,j]*(-2.0*u[i+1,j]+u[i,j]+8*(u[i+1,j]-u[i-1,j])+u[i-2,j])/(12.0*dxi) +
abs(u[i, j]) * (2.0*u[i+1,j]-u[i,j]-4.0*u[i+1,j]+6.0*u[i,j]-4.0*u[i-1,j]+u[i-2,j])/(4.0*dxi))
vrhs[i,j] -= exp(-(i-1)*dxi) * ( u[i,j]*(-2.0*v[i+1,j]+v[i,j]+8*(v[i+1,j]-v[i-1,j])+v[i-2,j])/(12.0*dxi) +
abs(u[i, j]) * (2.0*v[i+1,j]-v[i,j]-4.0*v[i+1,j]+6.0*v[i,j]-4.0*v[i-1,j]+v[i-2,j])/(4.0*dxi))
else
urhs[i,j] -= exp(-(i-1)*dxi) * ( u[i,j]*(-u[i+2,j]+8*(u[i+1,j]-u[i-1,j])+u[i-2,j])/(12.0*dxi) +
abs(u[i, j]) * (u[i+2,j]-4.0*u[i+1,j]+6.0*u[i,j]-4.0*u[i-1,j]+u[i-2,j])/(4.0*dxi))
vrhs[i,j] -= exp(-(i-1)*dxi) * ( u[i,j]*(-v[i+2,j]+8*(v[i+1,j]-v[i-1,j])+v[i-2,j])/(12.0*dxi) +
abs(u[i, j]) * (v[i+2,j]-4.0*v[i+1,j]+6.0*v[i,j]-4.0*v[i-1,j]+v[i-2,j])/(4.0*dxi))
end
end
end
#移流項 theta方向
for i=2:IM-1
for j=1:JM
if j==1
urhs[i,j] -= exp(-(i-1)*dxi) * ( v[i,j]*(-u[i,j+2]+8.0*(u[i,j+1]-u[i,JM-1])+u[i,JM-2])/(12.0*dtheta) +
abs(v[i,j]) * (u[i,j+2]-4.0*u[i,j+1]+6.0*u[i,j]-4.0*u[i,JM-1]+u[i,JM-2]) / (4.0*dtheta) )
vrhs[i,j] -= exp(-(i-1)*dxi) * ( v[i,j]*(-v[i,j+2]+8.0*(v[i,j+1]-v[i,JM-1])+v[i,JM-2])/(12.0*dtheta) +
abs(v[i,j]) * (v[i,j+2]-4.0*v[i,j+1]+6.0*v[i,j]-4.0*v[i,JM-1]+v[i,JM-2]) / (4.0*dtheta) )
elseif j==2
urhs[i,j] -= exp(-(i-1)*dxi) * ( v[i,j]*(-u[i,j+2]+8.0*(u[i,j+1]-u[i,j-1])+u[i,JM-1])/(12.0*dtheta) +
abs(v[i,j]) * (u[i,j+2]-4.0*u[i,j+1]+6.0*u[i,j]-4.0*u[i,j-1]+u[i,JM-1]) / (4.0*dtheta) )
vrhs[i,j] -= exp(-(i-1)*dxi) * ( v[i,j]*(-v[i,j+2]+8.0*(v[i,j+1]-v[i,j-1])+v[i,JM-1])/(12.0*dtheta) +
abs(v[i,j]) * (v[i,j+2]-4.0*v[i,j+1]+6.0*v[i,j]-4.0*v[i,j-1]+v[i,JM-1]) / (4.0*dtheta) )
elseif j==JM
urhs[i,j] -= exp(-(i-1)*dxi) * ( v[i,j]*(-u[i,3]+8.0*(u[i,2]-u[i,j-1])+u[i,j-2])/(12.0*dtheta) +
abs(v[i,j]) * (u[i,3]-4.0*u[i,2]+6.0*u[i,j]-4.0*u[i,j-1]+u[i,j-2]) / (4.0*dtheta) )
vrhs[i,j] -= exp(-(i-1)*dxi) * ( v[i,j]*(-v[i,3]+8.0*(v[i,2]-v[i,j-1])+v[i,j-2])/(12.0*dtheta) +
abs(v[i,j]) * (v[i,3]-4.0*v[i,2]+6.0*v[i,j]-4.0*v[i,j-1]+v[i,j-2]) / (4.0*dtheta) )
elseif j==JM-1
urhs[i,j] -= exp(-(i-1)*dxi) * ( v[i,j]*(-u[i,2]+8.0*(u[i,j+1]-u[i,j-1])+u[i,j-2])/(12.0*dtheta) +
abs(v[i,j]) * (u[i,2]-4.0*u[i,j+1]+6.0*u[i,j]-4.0*u[i,j-1]+u[i,j-2]) / (4.0*dtheta) )
vrhs[i,j] -= exp(-(i-1)*dxi) * ( v[i,j]*(-v[i,2]+8.0*(v[i,j+1]-v[i,j-1])+v[i,j-2])/(12.0*dtheta) +
abs(v[i,j]) * (v[i,2]-4.0*v[i,j+1]+6.0*v[i,j]-4.0*v[i,j-1]+v[i,j-2]) / (4.0*dtheta) )
else
urhs[i,j] -= exp(-(i-1)*dxi) * ( v[i,j]*(-u[i,j+2]+8.0*(u[i,j+1]-u[i,j-1])+u[i,j-2])/(12.0*dtheta) +
abs(v[i,j]) * (u[i,j+2]-4.0*u[i,j+1]+6.0*u[i,j]-4.0*u[i,j-1]+u[i,j-2]) / (4.0*dtheta) )
vrhs[i,j] -= exp(-(i-1)*dxi) * ( v[i,j]*(-v[i,j+2]+8.0*(v[i,j+1]-v[i,j-1])+v[i,j-2])/(12.0*dtheta) +
abs(v[i,j]) * (v[i,j+2]-4.0*v[i,j+1]+6.0*v[i,j]-4.0*v[i,j-1]+v[i,j-2]) / (4.0*dtheta) )
end
end
end
#時間発展
for i=2:IM-1
for j=1:JM
u[i, j] += dt * urhs[i, j]
v[i, j] += dt * vrhs[i, j]
end
end
end
VelocityEQ (generic function with 1 method)
const IM=81
const JM=121
const MAXITP=10000
const OMEGA=1.0
const ERRORP=1e-4
const MAXSTEP=10000
const RE=100
#const IO_FREQUENCY=50
100
function main(;step=10000,freq=100)
CFL=0.2
dxi=0.05
dtheta=2*pi/(JM-1)
dt=CFL*min(dxi,dtheta)
rho=1.0
p0=0.0
u=zeros(IM,JM)
v=zeros(IM,JM)
p=zeros(IM,JM)
rhs=zeros(IM,JM)
urhs=zeros(IM,JM)
vrhs=zeros(IM,JM)
vx=zeros(IM,JM)
vy=zeros(IM,JM)
zeta=zeros(IM,JM)
Ps=Array{Array{Float64,2}}(undef, step÷freq+1)
Us=Array{Array{Float64,2}}(undef, step÷freq+1)
Vs=Array{Array{Float64,2}}(undef, step÷freq+1)
Zs=Array{Array{Float64,2}}(undef, step÷freq+1)
#初期条件
for i=1:IM
for j=1:JM
u[i,j]=(1.0 - exp(-2*(i-1)*dxi))*cos((j-1)*dtheta)
v[i,j]=-(1.0 + exp(-2*(i-1)*dxi))*sin((j-1)*dtheta)
p[i,j]=p0-rho/2.0*(u[i,j]^2 + v[i,j]^2)
end
end
BCforP(p)
BCforV(u,v,dtheta)
num=1
Ps[num]=copy(p)
for j=1:JM
for i=1:IM
vx[i,j]=u[i,j]*cos((j-1)*dtheta) - v[i,j]*sin((j-1)*dtheta)
vy[i,j]=u[i,j]*sin((j-1)*dtheta) + v[i,j]*cos((j-1)*dtheta)
end
end
Us[num]=copy(vx)
Vs[num]=copy(vy)
for j=1:JM
for i=1:IM-1
vxp=u[i,ifelse(j==JM,2, j+1)]
vxm=u[i,ifelse(j==1,JM-1,j-1)]
dvx=(vxp-vxm)/(2*dtheta)
vt =v[i,j]
if i==1
dvt=(-3*v[1,j]+4*v[2,j]-v[3,j])/(2*dxi)
else
dvt=(v[i+1,j]-v[i-1,j])/(2*dxi)
end
zeta[i,j]=exp(-(i-1)*dxi)*(dvt+vt-dvx)
end
end
Zs[num]=copy(zeta)
# main loop
for itr=1:step
SolvePoisson(p, u,v, dtheta, dxi, dt,rhs)
BCforP(p)
VelocityEQ(p,u,v,dtheta,dxi,dt,urhs,vrhs)
BCforV(u,v,dtheta)
if itr%freq==0
num+=1
Ps[num] = copy(p)
for j=1:JM
for i=1:IM
vx[i,j]=u[i,j]*cos((j-1)*dtheta) - v[i,j]*sin((j-1)*dtheta)
vy[i,j]=u[i,j]*sin((j-1)*dtheta) + v[i,j]*cos((j-1)*dtheta)
end
end
Us[num]=copy(vx)
Vs[num]=copy(vy)
for j=1:JM
for i=1:IM-1
vxp=u[i,ifelse(j==JM,2, j+1)]
vxm=u[i,ifelse(j==1,JM-1,j-1)]
dvx=(vxp-vxm)/(2*dtheta)
vt =v[i,j]
if i==1
dvt=(-3*v[1,j]+4*v[2,j]-v[3,j])/(2*dxi)
else
dvt=(v[i+1,j]-v[i-1,j])/(2*dxi)
end
zeta[i,j]=exp(-(i-1)*dxi)*(dvt+vt-dvx)
end
end
Zs[num]=copy(zeta)
print(itr, " ")
end
end
return Ps,Us,Vs,Zs
end
main (generic function with 1 method)
@time Ps,Us,Vs,Zs=main(step=20000,freq=200);
200 400 600 800 1000 1200 1400 1600 1800 2000 2200 2400 2600 2800 3000 3200 3400 3600 3800 4000 4200 4400 4600 4800 5000 5200 5400 5600 5800 6000 6200 6400 6600 6800 7000 7200 7400 7600 7800 8000 8200 8400 8600 8800 9000 9200 9400 9600 9800 10000 10200 10400 10600 10800 11000 11200 11400 11600 11800 12000 12200 12400 12600 12800 13000 13200 13400 13600 13800 14000 14200 14400 14600 14800 15000 15200 15400 15600 15800 16000 16200 16400 16600 16800 17000 17200 17400 17600 17800 18000 18200 18400 18600 18800 19000 19200 19400 19600 19800 20000 32.324101 seconds (151.83 k allocations: 129.902 MiB, 0.08% gc time)
xi = range(0.0, 4.0, length=IM)
theta = range(0.0, 2*pi, length=JM)
r=exp.(xi)
xx = zeros(IM,JM); yy=zeros(IM,JM)
for i=1:IM
for j=1:JM
xx[i,j]=r[i]*cos(theta[j])
yy[i,j]=r[i]*sin(theta[j])
end
end
plev=range(-1.0, 0.75, length=31)
ulev=range(-0.5, 1.5, length=26)
vlev=range(-1, 1, length=21)
N=length(Ps)
anim = @animate for t=1:N
cfp=contourf(xx, yy, Ps[t],levels=plev,
xlim=(-5,15),
ylim=(-5,5),
clim=(-1,0.75),
color=:plasma,
aspect_ratio=1,
title="Pressure",
cbar=true,clabels=false
)
cfu=contourf(xx, yy, Us[t],levels=ulev,
xlim=(-5,15),
ylim=(-5,5),
clim=(-0.5,1.5),
color=:turbo,
aspect_ratio=1,
title="U",
cbar=true,clabels=false
)
cfv=contourf(xx, yy, Vs[t],levels=vlev,
xlim=(-5,15),
ylim=(-5,5),
clim=(-1,1),
color=:balance,
aspect_ratio=1,
title="V",
cbar=true,clabels=false
)
plot(cfp,cfu,cfv, size=(800,800*1.4),layout=(3,1))
end
filename="karman_enchu.gif"
gif(anim, filename, fps=10)
sys:1: UserWarning: The following kwargs were not used by contour: 'label' [ Info: Saved animation to C:\Users\****\Desktop\Julia\karman\karman_enchu.gif
zlev=range(-5,5,length=41)
contourf(xx, yy, Zs[101],levels=zlev,
xlim=(-5,15),
ylim=(-5,5),
clim=(-5,5),
color=:curl,
aspect_ratio=1,
title="Curl",
cbar=true,clabels=false
)
anim = @animate for t=1:N
cfz=contourf(xx, yy, Zs[t],levels=zlev,
xlim=(-5,15),
ylim=(-5,5),
clim=(-5,5),
color=:curl,
aspect_ratio=1,
title="Curl",
cbar=true,clabels=false
)
plot(cfz, size=(800,800*1.4/3))
end
filename="curl_enchu.gif"
gif(anim, filename, fps=10)
[ Info: Saved animation to C:\Users\****\Desktop\Julia\karman\curl_enchu.gif